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ABSTRACT 


Finite-difference,  time-domain  (FDTD)  techniques  hold  much  promise  for  performing  realistic  simu¬ 
lations  of  sound  propagation  through  complex,  dynamic  outdoor  environments.  This  report  focuses  on  a 
key  aspeet  of  FDTD  in  the  atmosphere,  namely  the  incorporation  of  a  moving  background  medium  (wind 
and  turbulence  in  the  atmosphere)  into  the  calculations.  Appropriate  differential  equations  for  FDTD 
simulation  of  sound  propagation  in  a  moving  fluid  are  discussed.  It  is  shown  that  FDTD  calculations  are 
not  possible  with  this  equation  set  when  using  the  staggered  grid,  “leap-frog”  approach,  which  is  com¬ 
mon  for  FDTD  simulation  of  other  types  of  wave  propagation.  Various  finite-difference  operators  that  are 
valid  for  a  moving  medium,  such  as  Runge-Kutta  and  an  unstaggered  leap-frog  approach,  are  discussed 
and  compared.  It  is  shown  that  a  rigorous  FDTD  solution  in  a  moving  medium  requires  storing  the  field 
variables  over  at  least  two  time  steps,  thereby  requiring  at  least  twice  as  much  computer  memory  as  the 
customary  staggered  grid.  Several  other  topies  pertinent  to  FDTD  simulation  of  sound  propagation  in  the 
atmosphere  are  discussed,  including  implementation  of  porous  ground  layers,  absorbing  boundaries,  and 
rigid  surfaces.  Example  calculations  demonstrate  the  performance  of  the  various  finite-difference  opera¬ 
tors  for  a  high  Mach  number,  uniform  flow.  Other  example  calculations  show  FDTD  calculations  for 
propagation  above  rigid  and  porous  ground  surfaces,  over  rigid  barriers,  and  through  turbulence.  With 
sufficiently  dense  spatial  grids,  very  good  agreement  can  be  obtained  between  the  FDTD  calculations  and 
known  theoretical  solutions. 
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All  product  names  and  trademarks  cited  are  the  property  of  their  respective  owners.  The  findings  of  this  report  are  not  to 
be  construed  as  an  official  Department  of  the  Army  position  unless  so  designated  by  other  authorized  documents. 
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1  INTRODUCTION 

Acoustical  sensing  systems  are  expected  to  play  a  key  role  in  the  Army’s  Future 
Force  by  providing  non-line-of-sight  surveillance  over  wide  areas.  It  is  well  known  that 
the  performance  of  battlefield  acoustical  sensors  is  affected  by  complex  sound  propaga¬ 
tion  phenomena  occurring  in  outdoor  settings,  such  as  reflections  from  trees  and  build¬ 
ings,  ground  interactions,  scattering  by  turbulence,  refraction  by  wind  and  temperature 
gradients,  and  diffraction  over  hills.  Comprehensive,  controlled  field  studies  of  the  effects 
of  these  phenomena  on  sensors  are  extremely  difficult  and  expensive  to  perform.  There¬ 
fore,  a  high-quality  simulation  capable  of  reproducing  these  phenomena  would  be 
immensely  valuable  for  the  development  of  effective,  reliable  acoustical  sensing  systems. 

Most  currently  used  numerical  methods  for  outdoor  sound  propagation,  such  as  the 
fast  field  program  (FFP)  and  parabolic  equation  (PE),  are  incapable  of  simulating  all  of 
the  propagation  phenomena  mentioned  above.  [The  reader  may  refer  to  Attenborough  et 
al.  (1995)  and  Salomons  (2001)  for  detailed  discussions  of  the  FFP  and  PE  methods.]  On 
the  other  hand,  many  of  these  phenomena  are  readily  handled  with  finite-difference,  time- 
domain  (FDTD)  techniques.  Reflections  from  and  diffraction  around  objects  such  as 
buildings,  trees,  and  hills  can  be  readily  implemented  with  the  FDTD  approach,  as  can 
dynamic  source  inputs  and  turbulent  scattering.  But  there  are  significant  drawbacks  of 
FDTD  that  have  so  far  prevented  its  widespread  use  in  outdoor  sound  propagation  calcu¬ 
lations.  In  particular,  it  is  very  computationally  intensive.  Other  inherent  difficulties 
include  implementing  a  time-domain  impedance  boundary  condition  for  the  ground  and 
controlling  numerical  instabilities  at  sudden  contrasts  in  material  properties.  But,  given 
the  inherent  flexibility  of  the  FDTD  technique  and  the  increasing  capability  of  modem 
computers,  FDTD  should  play  an  increasingly  prominent  role  in  outdoor  sound  propaga¬ 
tion  modeling. 

Recent  papers  by  Blumrich  and  Heimann  (2002)  and  Salomons  et  al.  (2002)  present 
2-D  FDTD  calculations  for  the  atmosphere  that  include  the  effects  of  motion  (wind  and 
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turbulence)  in  the  propagation  medium  and  address  issues  related  to  ground  interactions. 
Liu  and  Moran  (2002)  demonstrated  3-D  FDTD  simulations  of  propagation  in  a  non¬ 
moving  atmosphere  with  a  barrier  and  a  dense  stand  of  trees.  Because  of  the  substantial 
demands  on  computer  memory,  the  3-D  simulations  were  run  on  a  distributed,  parallel¬ 
processing  computer. 

In  this  report  we  discuss  in  detail  the  implementation  of  FDTD  for  sound  propagation 
in  a  moving  atmosphere.  By  “moving,”  we  mean  that  motion  in  the  medium  through 
which  the  sound  propagates  has  a  significant  impact  on  the  propagation,  which  is  gener¬ 
ally  true  when  the  typical  speed  of  motion  in  the  medium  is  not  very  much  smaller  than 
the  speed  of  wave  propagation  through  the  medium.  Motion  in  the  propagation  medium 
has  been  justifiably  neglected  in  previous  FDTD  simulations,  including  seismic  propaga¬ 
tion  in  the  earth  and  electromagnetic  propagation  in  the  atmosphere,  but  it  cannot  be 
neglected  for  sound  propagation  in  the  atmosphere.  For  example,  the  speed  of  light  in  air, 
about  3  X  10*  m  s“\  is  thirty  million  to  one  hundred  million  times  larger  than  typical  wind 
speeds  in  the  atmosphere  (3-10  m  s“').  The  speed  of  sound,  about  340  m  s“\  is  only 
about  34-1 10  times  larger.  As  a  result,  wind  and  turbulence  in  the  atmosphere  have 
important  effects  on  sound  propagation. 

In  Section  2  of  this  paper,  we  discuss  the  equations  for  sound  propagation  in  a  mov¬ 
ing  medium.  The  spatially  staggered  finite-difference  implementation  of  this  equation  set, 
in  which  the  acoustic  pressure  and  particle  velocities  are  stored  at  different  locations  in 
space,  is  described  in  Section  3.  Section  4  discusses  advancement  of  the  solution  in  time, 
on  grids  that  are  both  staggered  and  unstaggered  in  time.  We  point  out  that  the  conven¬ 
tional  staggered-in-time,  “leap-frog”  approach  for  a  nonmoving  medium  cannot  be 
applied  directly  to  a  moving  medium.  The  unstaggered-in-time  approach  allows  widely 
practiced  numerical  techniques,  such  as  Runge-Kutta  integration,  to  be  readily  applied  to 
the  problem.  In  Section  5,  several  other  topics  pertinent  to  FDTD  simulation  of  sound 
propagation  in  the  atmosphere  are  discussed,  including  implementation  of  porous  ground 
layers,  absorbing  boundaries,  and  rigid  surfaces.  Lastly,  some  example  calculations 
involving  propagation  through  a  moving  atmosphere  are  presented  in  Section  6. 
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2  EQUATIONS  FOR  SOUND  PROPAGATION 
IN  A  MOVING  MEDIUM 

Most  currently  used  techniques  for  calculating  sound  propagation  in  the  atmosphere, 
such  as  the  fast  field  program  and  parabolic  equation  (FFP  and  PE),  are  based  on  solution 
of  the  wave  equation  or  its  parabolic  approximation.  The  wave  equation  is  second  order 
with  respect  to  its  differential  operators  in  space  and  time.  FDTD,  however,  is  more  read¬ 
ily  applied  to  first-order  partial  differential  equations.  The  following  coupled,  first-order 
equations  for  the  acoustic  pressure  p  and  acoustic  particle  velocity  w  (Ostashev  1987, 
1997,  Aldridge  2002)  serve  as  an  appropriate  starting  point  for  FDTD  sound  propagation 
calculations  in  a  moving  medium: 

^  = -(v  V)j9-/7C^V-W-|-/7C^g, 

^  w  •  V)  v  -  ( V  •  V)w  -  — -t  F//7.  (2) 

Here,  p  is  ambient  medium  density,  c  is  the  adiabatic  speed  of  sound,  and  v  is  the 
wind  velocity.  The  quantities  F  and  Q  represent  sources:  the  former  is  a  force  acting  on 
the  medium,  whereas  the  latter  is  a  mass  source.  (These  may  be  interpreted  as  dipole  and 
monopole  pressure  sources,  respectively.)  As  is  typical  in  aeroacoustics,  these  equations 
assume  that  (1)  the  sound  wave  is  a  small  perturbation  to  the  background  state  of  the 
medium,  (2)  the  medium  behaves  adiabatically  (i.e.,  the  material  derivative  of  the  specific 
entropy  of  the  background  medium  vanishes),  and  (3)  the  divergence  of  the  background 
velocity  field  (the  turbulence)  is  zero.  The  gradient  of  the  background  pressure  is  also 
neglected;  this  approximation  is  reasonable  for  sound  propagation  near  the  ground  but 
would  need  to  be  reconsidered  for  inffasound  propagating  between  the  ground  and  upper 
atmosphere. 

The  terms  (v  ■  V  )p  and  ( v  •  V )  w  in  Eq.  1  and  2  represent  the  transport  of  the 
sound  wave  by  the  atmospheric  wind.  These  are  commonly  called  advective  terms.  The 
term  (  w  •  V  )  v  in  Eq.  2  represents  enhancement/diminishment  of  the  acoustic  particle 
velocity  by  spatial  variations  in  the  wind  field.  For  a  nonmoving  medium  (v  =  0),  one  can 
derive  the  customary  wave  equation  for  the  acoustic  pressure  by  taking  the  time  deriva¬ 
tive  of  Eq.  1  and  the  divergence  of  Eq.  2.  However,  the  additional  terms  involving  the 
wind  velocity  introduce  significant  complications  when  attempting  to  derive  wave  equa¬ 
tions  or  parabolic  equations  for  use  in  an  FFP  or  PE  (Pierce  1990,  Ostashev  et  al.  1997, 
Eingevitch  et  al.  2002).  An  FDTD  code  based  on  the  first-order  Eq.  1  and  2  would  be 
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more  general  and  accurate  than  most  current  sound  propagation  formulations,  despite  the 
fairly  simple  appearance  of  the  differential  equations. 

In  this  report  we  consider  the  FD  implementation  of  Eq.  1  and  2  in  two  dimensions. 
This  simplifies  the  presentation  considerably  and  makes  it  possible  to  run  calculations  on 
single-processor  computers.  The  generalization  of  the  results  to  three  dimensions  is  quite 
straightforward  but  more  computationally  intensive.  In  two  dimensions,  Eq.  1  and  2 
become 


dp 

8t 


(  d  d 
V^dx^^^dy 


P- 


^  dx 


+ 


dWy  '' 

dy  ^ 


(3) 


dt 


d  d^ 

Wy  —  +  W,,  —  V,- 

^dx  ^  dy) 


+  bF^, 


(4) 


and 


dWy 

dt 


d  d) 
^  dy 


f  d  d 
V^dx^^y  dy 


+  bFy, 


(5) 


where  b  =  Mp  is  the  mass  buoyancy  and  c  =  pc'  is  the  adiabatic  bulk  modulus. 
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3  STAGGERED  SPATIAL  FINITE-DIFFERENCE  GRID 

In  this  section  we  consider  second-order  approximations  of  the  spatial  derivatives  in 
Eq.  3-5.  The  spatial  finite-difference  analysis  in  this  report  is  based  entirely  on  a  grid  for 
the  pressure  and  particle  velocities  that  is  staggered  in  space,  as  shown  in  Figure  1 .  In  this 
grid  the  pressure  is  stored  at  whole  nodes,  namely  x  =  iAx  and  y  =j^y,  where  i  and  j  are 
integers  and  Ax  and  Ay  are  the  grid  spacings  in  the  x  and  y  directions.*  The  x  components 
of  the  acoustic  particle  velocity,  w„  are  staggered  (offset)  by  Ax/2  in  the  x  direction.  The 
y  components  of  the  acoustic  particle  velocity,  Wy,  are  staggered  by  Ay/2  in  the  y  direc¬ 
tion.  An  additional  convention  in  this  report  is  that  the  quantities  b,  k,  and  Q  are  stored  at 
the  pressure  nodes.  The  quantities  v^:  and  At  are  stored  at  the  w^:  nodes,  and  Vy  and  A,,  are 
stored  at  the  Wy  nodes. 


►  w, 

A 

•  P 


Figure  1.  Spatially  staggered  finite-difference  grid  for  solution 
of  acoustic  wave  propagation. 


Centered,  finite-difference  solution  of  Eq.  3-5  requires  an  evaluation  of  the  right- 
hand  sides  of  these  equations  at  the  grid  nodes  where  the  field  variable  is  stored.  For 


*The  spatially  staggered  grid  in  Fig.  1  matehes  the  one  normally  used  for  wave  propagation 
ealeulations  in  nonmoving  media.  Given  the  diffieulties  of  applying  the  staggered  grid  to  a 
moving  medium,  it  would  be  worthwhile  eonsidering  the  potential  advantages  of  spatially 
unstaggered  grids  in  this  eontext.  Sueh  eonsideration  is  not  given  in  this  report,  however. 
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example,  we  will  need  fmite-differenee  approximations,  centered  atx  =  iAx  and  j 
for  all  terms  on  the  right-hand  side  of  Eq.  3.  Finite  difference  approximations  centered  on 
X  =  (f  +  1/2) Ax  and  j  =J^y  are  needed  for  Eq.  4,  whereas  approximations  centered  onx  = 
iAx  and  j  =  (/  +  1/2) Aj  are  needed  for  Eq.  5. 

Consider  first  the  source  terms  and  coefficients  involving  the  medium  properties.  The 
source  terms  can  be  used  directly,  since  they  are  already  stored  at  the  grid  nodes  where 
the  ED  approximations  are  centered.  The  same  is  true  of  k,  which  is  stored  at  the  pressure 
grid  nodes  and  needed  in  Eq.  3.  Regarding  Eq.  4  and  5,  the  values  for  b  can  be  deter¬ 
mined  at  the  needed  locations  by  averaging  neighboring  values: 

+  b  [  iAx,  jAy,  t 

(6) 


b\_{i  +  \l  2)  Ax,  jAy,  t  ]  = 


b{{i  +  \)Ax,jAy,t] 


One  of  the  main  motivations  for  using  the  spatially  staggered  grid  is  that  it  conven¬ 
iently  provides  centered  spatial  differences  for  many  of  the  derivatives  in  Eq.  3-5.  In 
particular, 

dw,^  ( iAx,  jAy,  t)  ^  w^\_{i  +  \l  2)  Ax,  jAy,  -  w^\_{i  -  M  2)  Ax,  jAy,  t  ] 

&  "  aJc  ’  * 

dWy  ( iAx,  jAy, t  )  [/Ax,(  j  +  1/2)  Ay,t  ]  -  Wy  [/Ax,(  j  -  l/2)Ay,t  ] 

^  ""  A^^  ’  * 


dp\_{i  +  l/2)Ax,7Ay,t  ]  r[(/  +  l)Ax,7Ay,t]  -  p\_iAx,jAy,t'\ 

dx  Ax 


(10) 


and 


dp\_iAx,{j  +  1/2) Ay,t ]  p\_iAx,{j  +  1) Ay,t ]  -  p\_iAx,jAy,t'\  ^ 

dy  Ay 

At  this  point  we  have  developed  centered  approximations  for  all  of  the  terms  in  Eq. 
3-5  that  are  present  for  a  nonmoving  medium.  Implementation  of  the  remaining  terms  is 
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somewhat  more  complieated,  however.  For  example,  the  derivatives  of  the  pressure  field 
in  Eq.  3,  dpidx  and  Bp! By,  cannot  be  centered  atx  =  /Ax  and  j  from  approximations 
across  a  single  grid  interval.  We  can  develop  centered  approximations  across  two  grid 
intervals,  however: 


Bp{iAx,jAy,t  )  ^  p\_{i  +  ^)Ax,jAy,t'\  -/>[(/-  \  )Ax,jAy,t'\ 
Bx  2Ax 


and 


Bp  ( /Ax,  jAy,  t ) 
By 


p[iAx,{j  +  i)Ay,t]- p[iAx,{j -l)Ay,t] 
2Ay 


(12) 


(13) 


To  find  the  wind  velocity  components  and  Vy  evaluated  atx  =  /Ax  andy  =7 Ay 
(multiplying  the  derivatives  BpIBx  and  BpIBy,  respectively,  in  Eq.  3),  we  average  values  at 
neighboring  grid  points: 


v^{iAx,jAy,t) 


[(  /  +  1/2) Ax,7Ay,/]  +  v^\_{i  - 1/2) Ax,7Ay,/] 
2 


(14) 


v^[Ax,(7  +  l/2)Ay,/]  +  v^[/Ax,(7-l/2)Ay,/] 
Vy{iAx,jAy,t)=^ - - - 

A  similar  approach  to  Eq.  4  and  5  leads  to  the  approximations 

[  ( /  +  1  /  2 )  Ax,  7  Ay,  t  ] 

Bx 


[  ( /  +  3  /  2 )  Ax,  7  Ay,  t^  -  -  M  2)  Ax,  7  Ay,  t  ] 

2Ax 


(15) 


(16) 


Bw^  [  ( /  +  1  /  2 )  Ax,  7  Ay,  t  ] 

By 

[(/  +  l/2)Ax,(7  +  l)Ay,/]- [(/  +  l/2)Ax,(7-l)Ay,/] 

2Ay 


w. 


(17) 
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8Wy  [  /Ax,  (  7  +  1  /  2  )  Aj,  t  ] 

& 

[(  /  +  1)Ax,(7  +  l/2)Aj,/]  -  w^[(  /  -  1)Ax,(7  +  l/2)Aj,/] 

2Ax 


(18) 


and 


[  /Ax,  (  7  +  1  /  2  )  Aj,  /  ] 

dy 

Wy  [  /Ax,  (  7  +  3  /  2  )  Aj,  /  ]  -  Wy  [  /Ax,  (  7  -  1  /  2  )  Ay,  t  ] 
“  2A^^ 


(19) 


The  same  approximations  hold  for  derivatives  of  the  wind  field  when  Wx  is  replaced 
by  Vx  and  Wy  by  Vy.  In  Eq.  4  the  quantities  Wy  and  Vy  (multiplying  the  derivatives  dvjdy 
and  dwjdy,  respectively)  are  needed  at  the  grid  point  x  =  (/  +  1/2) Ax  andy  =7 Ay.  Refer¬ 
ring  to  Figure  1 ,  a  reasonable  way  to  obtain  these  quantities  would  be  to  average  the  four 
closest  grid  nodes: 


[  ( /  +  1  /  2  )  Ax,  7 Ay,  t  ] 

==  +  1)'^>(7  +  l/2)Ay,/]  +  Wy  [iAx,(j  +  l/2)Ay,/]  (20) 

+Wy  [(/  +  1)Ax,(7  -  l/2)Ay,/]  +  Wy  [iAx,(j  -  l/2)Ay,/]j- , 

and  likewise  for  Vy.  Similarly,  in  Eq.  5  the  quantities  Wx  and  Vx  (multiplying  the  deriva¬ 
tives  dvyidx  and  dwyldx,  respectively)  are  needed  at  the  grid  point  x  =  /Ax  and 
y  =  (/  +  1/2) Ay.  We  have 

Wx  [  /Ax,  (  7  +  1  /  2  )  Ay,  /  ] 

~  4  [( ^  ■*“  2)  Ax,  (  7  +  1 )  Ay,  /  ]  +  Wx\_{i  +  1/2  )  Ax,  7  Ay,  /  ]  (21) 

[  ( ^  “  1  /  2  )  Ax,  (  7  +  1 )  Ay,  /  ]  +  w^[(/-l/2)  Ax,  7  Ay,  /  ]  , 


and  likewise  for  v^. 
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Equations  6-2 1  provide  finite-difference  approximations  for  all  of  the  quantities 
appearing  on  the  right-hand  sides  of  Eq.  3-5.  We  next  turn  our  attention  to  finite- 
difference  approximations  for  the  temporal  derivatives. 


10 


ERDC/CRRELTR-04-12 


4  TEMPORAL  FINITE-DIFFERENCE  GRIDS 

In  this  section  we  consider  two  alternative  approaches  to  advancing  the  wavefield 
variables  in  time.  The  first  is  based  on  an  unstaggered-in-time  grid,  meaning  that  the 
pressure  and  acoustic  particle  velocities  are  stored  on  the  same  time  level.  The  second  is 
based  on  a  staggered-in-time  grid,  in  which  the  pressure  and  particle  velocities  are  stored 
on  alternating  half  time  steps. 


Solution  on  a  Grid  Unstaggered  in  Time 

Let  us  define  fp  as  the  right-hand  side  of  Eq.  3.  We  write 


dp  ( i/Sx,  jky,  t ) 
dt 


=  4[«Ax,7Aj,p(0,w^(0,w^(0,s(0], 


(22) 


where  p(t),  Wx(05  and  w/t)  represent  matrices  containing  the  pressures  and  particle 
velocities  at  all  available  grid  nodes.  Furthermore,  s(t)  represents  the  source  and  medium 
properties  {b,  k,  v„  Vy,  Q,  F„  and  Fy)  at  all  available  grid  nodes.*  From  Eq.  8-9  and  12- 
15,  the  second-order  FD  approximation  of fp  is 

/p[/Ar,7Aj,p(0,w^(0,w^(t),s(0]  = 

1 


4Ax 


{v^[{i  +  \/2)Ax,jAy,t]  +  v^[{i-\/2)Ax,jAy,t]} 


4Aj 


X  { p  [  ( /  -I- 1 )  Ax,  jAy,  t]-  p[{i  -1)  Ax,  JAy,  t  ] } 


[vy[Ax,{j  +  \l2)Ay,t]  +  Vy[iAx,{j -\l2)Ay,t]] 


(23) 


^{p[i^x,{j  +  \)Ay,t]-  p[iAx,{j  -\)Ay,t]} 

w^\_{i  +  \l  2)  Ax,  jAy,  t'\-w^\_{i-\l2)  Ax,  JAy,  t  ] 


-K{iAx,jAy,t) 


{■ 


Ax 


-!-■ 


Wy  [  iAx,  {j  +  \  1 2)  Ay,  t'\-Wy[  iAx,  (J  -1/2)  Ay,  t  ] 


At 


-r/r  ( iAx,  jAy,  t)Q(  iAx,  jAy,  t ) . 


} 


*Note  that ^[iAx,^^^,  p(t),  wft),  Wy(i),  s(t)]  in  actuality  depends  only  on  the  fields  at  grid  points 
in  the  immediate  vicinity  of  (;Ax,  jAy)  when  second-order  spatial  differencing  is  used.  The 
notation  used  here,  which  is  convenient  for  discussion  purposes,  is  general  enough  to  accommodate 
spatial  differencing  of  an  arbitrarily  high  order  as  well  as  staggered  and  unstaggered  spatial  grids. 
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Similarly  to  Eq.  22  we  define  and  fy  as  the  right-hand  sides  of  Eq.  4  and  5: 


and 


dwy  [  /Ax,  ( 7  +  1  /  2 )  Ay,  t  ] 
dt 


[/Ax,( 7  +  1/2)  Ay,p(/),w^(/),w^(/),s(/)],  (25) 


where,  from  the  equations  in  Seetion  3, 

fx  [(/  +  1  /  2)  Ax,  7 Ay,  p  (/) ,  W;,  (/) ,  (/) ,  s  (/)]  - 

ny  [(/  +  1/2)  Ax,  7  Ay,  /] 


2Ax 


{Vx  [(/  +  3  /  2)  Ax,  7 Ay,  /]  -  y,  [(/  - 1  /  2)  Ax,  7 Ay,  /]} 


-  {wy  [(/  +  1)  Ax,  (7  +  1  /  2)  Ay,  /]  +  Wy  [/Ax,  (7  +  1  /  2)  Ay,  /] 

+Wy  [(/  +  1)  Ax,  (7  -  1  /  2)  Ay,  /]  +  Wy  [/Ax,  (7  -  1  /  2)  Ay,  /]} 

X  {V;,  [(/  +  1  /  2)  Ax,  (7  +  1)  Ay,  /]  -  y,  [(/  +  1  /  2)  Ax,  (7  -  1)  Ay,  /]} 
[(/  +  1/2)  Ax,  7 Ay,  /] 


2Ax 


-  {ly-  [(/  +  3/2)  Ax,  7 Ay,  /]  -  ly.  [(/  -1/2)  Ax,  7 Ay,  /]} 


-  ^  {v^;  [(/  +  1)  Ax,  (7  +  1  /  2)  Ay,  /]  +  Vy  [/Ax,  (7  +  1  /  2)  Ay,  /] 

+Vy  [(/  +  1)  Ax,  (7  - 1  /  2)  Ay,  /]  +  Vy  [/Ax,  (7  -  1  /  2)  Ay,  /]} 

X  {ly,  [(/  +  1  /  2)  Ax,  (7  +  1)  Ay,  /]  -  ny  [(/  +  1  /  2)  Ax,  (7  -  1)  Ay,  /]} 

1 

{b  [(/  +  1)  Ax,  7  Ay,  t\  +  b  [/Ax,  7  Ay,  /]} 


(26) 


2Ax 


X  {7?  [(/  +  1)  Ax,  7 Ay,  t\-  p  [/Ax,  7  Ay,  /]} 
b  [(/  +  1)  Ax,  7  Ay,  t\  +  b  [/Ax,  7  Ay,  /] 


E;[(/  +  l/2)Ax,7Ay,/] 


and 
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fy  [i^,  (7  + 1  /  2)  Aj,  p  (t) ,  (?) ,  (?) ,  s  (?)]  - 

-  ^  {wy  [(/  +  1  /  2)  Aj,  {j  +  1)  Aj,  ?]  +  Mi  [(?  +  1  /  2)  Aj,  jAy,  t] 

[(^  -1/2)  Ax,  (7  +  1)  Ay,  ?]  +  Mi  [(?  -1/2)  Ax,  jAy,  ?]} 

X  {v^  [(?  +  1)  Ax,  (7  +  1  /  2)  Ay,  t]  -  Vy  [(?  -  1)  Ax,  (7  +  1  /  2)  Ay,  ?]} 

— — - ’^^2 Ay — ^  ^  (2  +  2 / 2)  Aj, ?]  -  Vy  [?Ax, (7  -  1  / 2)  Ay, ?]} 

-  ^  {v;c  [(/  +  1  /  2)  Ax,  (7  +  1)  Aj,  ?]  +  V;,  [(?  +  1  /  2)  Ax,  JAy,  ?] 

+v^  [(?  - 1  /  2)  Ax,  (7  +  1)  Aj,  ?]  +  V;,  [(?  - 1  /  2)  Ax,  JAy,  ?]} 

X  |Mi  [(?  +  1)  Ax,  (7  +  1  /  2)  Ay,  ?]  -  Mi  [(?  - 1)  Ax,  (7  +  1/2)  Ay,  ?]} 

— — - ’  ^^2 Ay — ^  ^  (2  +  2  /  2)  Aj,  ?]  -  Mi  [?Ax,  (7  -  1  /  2)  Aj,  ?]} 

-  {b  [iAx,  (7  +  1)  Ay,  t\  +  b  [iAx,  jAy,  ?]} 

X  {7?  [?Ax,  (7  +  1)  Aj,  ?]  -  7?  [??ix,  JAy,  ?]} 

[,Ar,0  +  1/2)A;,,,]^ 


Suppose  now  that  p,  w^,  w^,,  and  s  are  stored  at  discrete  time  levels  ?  =  lAt.  A  simple 
approximation  to  the  time  derivatives  that  we  could  consider  is  the  forward-difference 
approximation,  namely, 

dp{iAx,JAy,lAt)  7?[?Ax, 7'Ai, (/ +  1 ) A? ]  -  73[?Ax,7'Ai, /A?] 
dl  “  A,  ’ 

and  similarly  for  dw^  [(?  +  1/2)  Ax,  JAy,  /A?]/  dt  and  dwy  [?Ax,  (7  +  1/2) Ay,  /A?]/  dt.  Solving 
for  the  variables  at  the  time  ?  =  (/  +  1)A?  gives 
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p  [  /Ax,  j/S.y,  (/  +  l)A/]  =  p[  /Ax,  7  Aj,  /A/  ] 

+A^p  [/Ax,7Aj,p(/A/),Wj^  (/A/),w^  (/A/),s(/A/)], 


(29) 


[ ( /  +  1  / 2 )  Ax, 7 Aj, (/  +  l)A/]  =  i4';f[(/  +  l/  2) Ax, yAj, /A/ ] 

+A//).  [(  /  +  1/2)Ax,7Aj,p(/A/),Wj^  (/A/),w^  (/A/),s(/A/)], 

and 


(30) 


Wy  [  /Ax,  ( 7  +  1  /  2 )  Ay,  ( /  +  1 )  A/  ]  =  [  /Ax,  ( 7  +  1  /  2 )  Ay,  I  At  ] 

+Atfy  [  /Ax,  ( 7  +  1  /  2 )  Aj,  p  ( /A/ ) ,  Wj-  ( /A/ ) ,  ( /A/ ) ,  s  ( /A/ )  ] . 


(31) 


This  method  of  advancing  the  variables  in  time  based  on  forward  differences  is  often 
called  Euler’s  method  (e.g..  Burden  et  al.  1981,  Kreyszig  1988).  Unfortunately,  as  dis¬ 
cussed  in  Kasahara  (1977),  this  method  is  unconditionally  unstable.  Even  when  the  time 
step  is  made  very  small,  the  error  will  grow  exponentially  and  eventually  become 
unacceptable. 

A  better  approach  is  to  center  the  approximations  for  the  time  derivatives  by  using 
fields  from  two  previous  time  steps: 


dp  ( /Ax,  7  Ay,  I  At )  p  [  /Ax,  7  Ay,  (/  +  l)A/]-7?[  /Ax,  7  Ay,  ( /  -  1 )  A/  ] 

~dt  "  TAt  ■  * 

Note  that  the  derivative  at  /  =  I  At  is  being  estimated  from  the  fields  at  /  =  (/  -  1)A/ 
and  /  =  (/  +  1)A/.  This  leads  to  the  updating  equation 


p [ /Ax, 7 Ay, (/  +  1)A/]  =  77[ /Ax, 7 Ay, ( /  - 1 )  A/ ] 

r  n 

+2Atfp  [/Ax,7Ay,p(/A/),Wjt  (/A/),w^  (/A/),s(/A/)J, 

and  similarly  for  and  Wy.  Kasahara  (1977)  calls  this  method  the  “leap-frog”  scheme. 
Unlike  the  Euler  scheme,  it  is  generally  stable.  The  leap-frog  scheme  was  one  of  the  most 
common  in  atmospheric  modeling  at  the  time  of  Kasahara’ s  writing.  As  will  be  shown 
later  with  examples,  it  is  also  a  satisfactory  scheme  for  calculating  sound  propagation  in  a 
moving  medium.  However,  the  field  variables  must  be  stored  over  two  time  steps, 
thereby  doubling  memory  requirements.  In  recent  decades  the  leap-frog  scheme  has  been 
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replaced  by  somewhat  more  complicated  ones  that  usually  provide  more  accurate  calcu¬ 
lations  with  little  additional  computational  burden. 

One  important  class  of  numerical  methods  that  is  in  widespread  usage  today  is  the 
Runge-Kutta  methods  (Burden  et  al.  1981,  Kreyszig  1988).  The  second-order  Runge- 
Kutta  method  (also  called  Heun’s  method  or  the  Improved  Euler  method)  can  be  moti¬ 
vated  in  an  intuitive,  heuristic  manner.  It  begins  in  the  same  way  as  Euler’s  method,  by 
forming  an  estimate  for  the  field  at  the  next  time  step  based  on  Eq.  29-3 1 .  These 
approximations  are  indicated  here  with  superscripted  asterisks,  i.e.. 


p*  [  /Ax,  7  Ay,  (/-tl)A/]  =  p[  iAx,  jAy,  I  At  ] 

+Atfp  \_iAx,jAy,’p{lAt),yi^  (/A/ ),w^  (/A/),s(/A/)] , 


(34) 


w*  [  ( /  -t  1  /  2 )  Ax,  jAy,  (/-tl)A/]  =  iV;f[(/-tl/2)  Ax,  JAy,  I  At  ] 
-tA//'^  [(/-t  l/2)Ax,7Ay,p(/A/),Wj£  (/A/),w^  (/A/ ),s(/A/ )], 


(35) 


w*y  [  /Ax,  ( 7  -t  1  /  2 )  Ay,  ( /  -t  1 )  A/  ]  =  [  /Ax,  ( 7  -t  1  /  2 )  Ay,  I  At  ] 

(36) 

+tfy  [/Ax,(7-tl/2)Ay,p(/A/),Wjc  (/A/),w^  (/A/),s(/A/)]. 

Next  p*[(/  -I- 1)  At] ,  w*[(/  +  1)  At],  and  w*  [(/  -I- 1)  At]  are  used  to  estimate  a 
value  for  the  derivative  functions  fp,fx,  and  fy  at  the  end  of  the  time  step,  e.g., 

fp  [/Ax,7Ay,p[(/-tl)A/],w^  [(/  -t  1)  A/],w^  [(/  -t  1)  A/],s[(/  -t  1)A/]]  = 

r  n 

fp  [iAx,7Ay,p*  [(/-tl)A/],w^(/-tl)A/],w*  [(/-t  l)A/],s[(/-tl)A/]]. 

Then,  in  updating  the  fields,  the  derivatives  at  the  beginning  and  end  of  the  time  step 
are  averaged  to  approximate  the  derivative  centered  on  the  time  step.  Specifically  the 
pressure  field  is  advanced  by  the  equation 

p  [  /Ax,  jAy,  (/-tl)A/]  =  7?[  /Ax,  jAy,  I  At  ] 

+  y{/p  [/Ax,7Ay,p(/A/),w^(/A/),W3,  (/A/),s(/A/)] 

+fp  [/Ax,7Ay,p*[(/-tl)A/],w*  [(/  +  l)A/],w*  [(/-t  1)  A/],s[(/-t  1)A/]]| . 


(38) 
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The  particle  velocities  are  advanced  in  the  same  way. 

The  fourth-order  Runge-Kutta  method  (Burden  et  al.  1981,  Kreyszig  1988)  is  also 
widely  used.  We  do  not  attempt  a  heuristic  explanation  of  it  here  but  rather  just  state  the 
equations  as  they  appear  in  the  notation  of  this  report: 

{iAx,jAy,lAt)  -  Atfp  \_iAx,jAy,i^{lAt),yfx  (/At ),w^  (/At ),s(/At )],  (39) 


( iAx,  jAy,  I  At )  -  Atfp  [/Ax,  jAy,  p  ( /At )  -t  p'  ( /At )  /  2, 

(/At)  -t  (/At  )/2,Wj,  (/At )  -t  (/At  )/2,s[(  /  -t  1/2)  At]], 


(40) 


p^  ( iAx,  jAy, I  At )  -  Atfp  [/Ax,  jAy,  p(/At)-tp^(/At)/2, 

W;c  (/At )  -t  (/At  )/2,w^  (/At)  -t  (/At)/2,s[(/  -t  1/2)  At]], 


(41) 


p^  ( iAx,  jAy,  I  At )  -  Atfp  [/Ax,  jAy,  p(/At)-tp^(/At), 
W;f  ( /At )  -t  ( /At ) ,  ( /At )  -t  ( /At ) ,  s  [  /At  ]], 

p  [  iAx,  jAy,  (/-tl)At]  =  j3[  iAx,  jAy,  I  At  ] 

1  I  />*  ( iAx,  jAy,  lAt)  +  2p^  { iAx,  jAy,  I  At)  1 
^\+2p^  { iAx,  jAy, I  At )  +  p^  { iAx,  jAy,  I  At )  j 


(42) 


(43) 


Equations  for  are  the  same  except  with  f  replacing yj,  and  (/  +  1/2) Ax  replacing 
iAx;  equations  for  Wy  are  the  same  except  with  f,  replacing yj,  and  (/  +  1/2)  Ay  replacing 

j^y- 


Solution  on  a  Grid  Staggered  in  Time 

A  common  approach  to  FDTD  solution  for  wave  propagation  in  a  nonmoving 
medium  is  to  use  a  staggered,  centered  finite-differences  grid  in  time.  This  approach  is 
discussed,  for  example,  by  Yee  (1966)  and  Botteldooren  (1994).  The  main  motivation  for 
this  approach  is  that  it  produces  a  solution  with  second-order  accuracy  with  regard  to  the 
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temporal  grid  resolution  At,  while  requiring  storage  of  the  wavefield  variables  at  only  a 
single  time  step. 

Consider  Eq.  3-5,  as  specialized  for  a  nonmoving  medium: 


dp 

dt 


(  dw^  dwy  ^ 

^  dx  dy  ^ 


+  kQ, 


(44) 


dt 


+  bF^ 


and 


dWy 

dt 


+  bFy. 


(45) 


(46) 


Suppose  now  that  the  pressure  is  stored  on  the  whole  time  steps  and  the  velocities  at 
the  half  time  steps.  By  forming  centered  FD  approximations  in  time,  we  have,  for  the 
pressure 


0/>[/A\:,yAj,(  /  +  l/2)At] 
dt 


-  /p[/Aic,yAj,W;c[(/  +  l/2)At],w^[(/  +  l/2)At],s[(/  +  l/2)At]] 

^  p[ilXxJlXy,[  I  +  l)At]  -  p[ilXx,jlXy,llXt\ 

At 


(47) 


Solving  for p[iAx,jAy,  (/  +  l)At],  we  have 


p  [  iAx,  jAy,  (/  +  l)At]  =  j3[  iAx,  jAy,  I  At  ] 

+  Atfp\_iAx,jAy,yi,;[{l  +  l/2)At],w^[(/  +  l/2)At],s[(/  +  l/2)At]]. 


(48) 


A  similar  process  for  the  velocities  yields 


[(  /  +  1/2) Ax,7Ay,( /  +  1/2) At ]  =  [(  /  +  l/2)Ax,7Ay,( /  - 1/2) At ] 

+  At/'^[(t  +  l/2)Ax,7Ay,p(/At),s(/At)], 


w. 


(49) 
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and 


[/Ax,(  j  +  1/2) Aj,( /  +  l/2)Af]  =  [/Ax,(  j  +  1/2) Aj,(/  -  1/2) Af  ] 

+  A^^  [/Ax,(  j  +  l/2)Aj,p(/Af),s(/Af )]. 


(50) 


Note,  for  this  scheme,  that  fp  is  needed  in  Eq.  48  only  on  the  half  time  steps,  and  that 
only  the  particle  velocities  (not  the  pressure)  are  needed  for  its  calculation.  Furthermore, 
in  Eq.  49  and  50,  fx  and  fy  are  needed  only  on  the  whole  time  steps  and  require  only  the 
pressure.  Therefore,  the  solution  can  be  advanced  in  time  by  solving  for  the  pressure  on  a 
whole  time  step,  solving  for  the  particle  velocities  on  the  next  half  time  step,  solving  for 
the  pressure  at  the  next  whole  time  step,  etc.  The  wavefield  variables  can  be  overwritten 
in  place  as  they  are  updated.  Because  the  centered  differences  are  accurate  to  second- 
order,  this  scheme  provides  good  accuracy  as  well  as  efficiency  in  its  use  of  computer 
memory.  This  staggered  scheme  is  usually  called  the  leap-frog  method,  which  unfortu¬ 
nately  is  the  same  term  used  by  Kasahara  (1977)  in  connection  with  the  unstaggered, 
centered,  second-order  scheme  (Eq.  33). 

Extending  consideration  now  to  the  moving  medium,  we  see  in  comparing  Eq.  3  and 
44  that  pressure  is  no  longer  absent  from  fp.  Similarly  Wx  and  Wy  are  required  in  the 
calculation  of fx  and  fy. 


p  [  /Ax,  7  Ay,  (/  +  1)A/]  =  j3[  /Ax,  y  Ay,  /A/  ] 


+  ^tfp 


/Ax,7Ay,p[(/  +  l/2)A/],Wjt[(/  +  l/2)A/], 
w^[(/  +  l/2)A/],s[(/  +  l/2)A/] 


(51) 


A  similar  process  for  the  velocities  yields 


Wx\_{i  +  1/2) Ax,7Ay,( /  +  1/2) A/]  =  Wx\_{i  +  l/2)Ax,7Ay,( /  - 1/2) A/] 
+  A/4  [(  /  +  l/2)Ax,7Ay,p(/A/),Wj^(/A/),w^  (/A/),s(/A/)], 


and 


Wy  [  /Ax,  ( 7  +  1  /  2 )  Ay,  (/  +  l/  2)A/]  =  w^[  /Ax,  ( 7  +  1  /  2 )  Ay,  (/-1/2)A/] 
+  A/4  [/Ax,(  7  +  l/2)Ay,p(/A/),Wjt  (/A/),w^  (/A/),s(/A/)]. 


(53) 
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As  a  result,  it  is  no  longer  possible  to  implement  the  basic  staggered-in-time  scheme 
considered  in  the  previous  paragraph.  One  work-around  that  might  be  considered  is  to 
substitute  {p[(/  +  l)At]  -i-  p[/At]}/2  where  p[(/  +  l/2)At]  is  needed  in  Eq.  48.  However, 
one  is  then  left  with  an  equation  for />[/Ax,7Aj,  (/  -l-  l)At]  that  cannot  be  solved  without 
simultaneously  knowing  the  pressure  at  all  locations  at  t  =  (/  -l-  l)At.  Therefore,  this  pro¬ 
cedure  does  not  provide  an  explicit  method  for  updating  the  pressure  field  at  each  time 
step. 

Another  possibility  would  be  to  simply  use  the  pressure  field  p(/At)  in  place  of  p[(/  + 
l/2)At  when  evaluating^,  and  Wx[(/  -  l/2)At]  and  yfy[{l  -  l/2)At]  in  place  of  and 

Wy(lAt)  when  evaluating and/v.  This  nonrigorous  procedure  is  somewhat  similar  to 
using  the  Euler  (forward-difference)  method  from  the  previous  section  to  evaluate  the 
advective  terms  while  maintaining  the  staggered-in-time  grid  approach  for  the  terms 
related  to  the  nonmoving  medium.  From  a  programming  standpoint,  one  simply  updates 
(overwrites)  p  based  on  calculation  of fp  from  the  currently  stored  p,  w^,  and  w_>,,  and  then 
updates  W;^  and  w,,  based  on  calculation  of and  fy  from  the  currently  stored  values  for  p, 
Wx,  and  Wy.  This  alternating  process  is  repeated  until  the  calculation  reaches  the  desired 
end.  The  calculations  in  Salomons  et  al.  (2002)  and  Blumrich  and  Heimann  (2002) 
appear  to  be  based  on  such  a  procedure.  Salomons  et  al.  mention  (p.  486)  that  they  first 
solve  for  the  particle  velocities  and  then  use  these  updated  velocities  to  calculate  the 
pressure,  cautioning  that  “Use  of  the  non-updated  velocities  gives  an  unstable  solution.” 
Thus,  the  authors  appear  to  be  using  a  temporally  staggered  implementation,  although  it 
is  not  entirely  clear  from  the  equations  presented  in  the  text. 

A  rigorous  scheme  for  performing  centered,  finite-difference  calculations  on  a 
staggered-in-time  grid  has  been  proposed  by  D.  Aldridge.*  This  scheme  involves  storing 
the  field  variables  back  two  steps  in  time,  rather  than  just  one,  as  has  been  considered 
thus  far.  Forming  a  centered  finite  difference  across  two  time  steps,  one  has,  instead  of 
Eq.  51, 

p\_iAx,jAy,{l  +  l)At  ] 

-  p\_iAx,jAy,{^l -\)At'\  (54) 

-t2 At/),  [  / Ax,  7 Ay,  p  ( / At ) ,  Wx  ( / At ) ,  ( / At ) ,  s  ( /At )  ] . 

The  problem  now  is  that  the  velocities  are  unknown  on  the  time  step  t  =  /At.  This  can 
be  overcome  by  averaging  the  velocities  from  the  neighboring  two  time  steps;  that  is. 


Personal  communication,  2003. 
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w^[(/  +  l/  2)Af]  +  w^[(/-l/2)Af] 
2 


and 


w^[(^  +  l/2)Af]  +  w^[(/-l/2)Af] 

)  2 

Similarly,  the  velocities  can  be  updated  with  the  equations 


(55) 


(56) 


w^\_{i  +  1/2) Aj,7Ay,( /  +  1/2) At ]  =  w^\_{i  +  M 2)/Ax,jAy,{l  - 3/2) At ] 

(t  +  l/2)Av,7Ay,p[(/  -  l/2)At],W;f  [(/  -  l/2)At], 
w^.[(/-l/2)At],s[(/-l/2)At] 


(57) 


Wy  [  iAx,  (J  +  1/2)  Ay,  (/  +  l/  2)At]  =  Wy[  iAx,  (J  +  1/2)  Ay,  (/-3/2)At] 


+2Atfy 


iAx,{  j  +  l/2)Ay,p[(  /  -  l/2)At], 
W;,[(/-l/2)At],w^[(/-l/2)At],s[(/-l/2)At] 


(58) 


and 


p[(/-l/2)At] 


p[/At]  +  p[(/-l)At] 
2 


(59) 


This  scheme  is  similar  to  the  unstaggered  leap-frog  method  discussed  in  the  previous 
section.  The  main  differences  are  (1)  the  averaging  of  the  field  variables  from  the 
neighboring  time  steps,  displaced  ±1/2  time  of  the  centered  approximation,  and  (2)  the 
alternating  update  of  the  pressure  and  particle  velocities. 
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5  IMPLEMENTATION  OF  REFLECTIVE  SURFACES 

In  this  section  we  consider  boundary  conditions  for  the  finite-difference  solution  as 
well  as  the  incorporation  of  a  physically  realistic  ground  layer. 

Rigid  Surfaces 

Hard  surfaces  such  as  buildings,  walls,  and  ice  can  often  be  modeled  as  acoustically 
rigid.  By  definition,  the  acoustic  particle  velocity  vanishes  in  the  direction  normal  to  the 
rigid  surface.  A  practical  method  for  handling  the  rigid  boundary  is  the  method  of 
images.  In  this  method  the  actual  medium  below  the  rigid  boundary  is  replaced  by  the 
mirror  reflection  of  the  medium  above  the  boundary.  An  image  source  is  placed  at  a 
distance  below  the  boundary  that  is  equal  to  the  distance  above  the  boundary  for  the 
actual  source.  Therefore,  the  solutions  of  the  differential  equations  above  and  below  the 
boundary  are  the  same,  except  that  they  are  reversed  in  the  coordinate  direction  normal  to 
the  boundary.  As  a  result,  the  acoustic  particle  velocities  in  the  two  solutions  exactly 
cancel  at  the  boundary.  Because  the  solution  above  the  boundary  satisfies  the  original 
differential  equation  and  the  desired  boundary  condition,  we  are  guaranteed  that  it  is  the 
correct,  unique  solution.  (The  solution  below  the  boundary  is  incorrect,  but  that  is  not  a 
concern.) 

Let  us  consider  the  finite-difference  equations  23-27  in  the  presence  of  a  rigid 
surface.*  The  surface  is  assumed  to  be  positioned  aty  =  (/  -  1/2) Ay  in  Figure  1.  In  the 
pressure  derivative  equation,  Eq.  23,  we  of  course  would  set  Wj,[/A>c,  (/  -  l/2)Ay,  t]  =  0. 
One  also  must  have  v^[/Aj,  (/  -  1/2) Ay,  t]  =  0  at  the  rigid  boundary.  (Otherwise,  the  wind 
fields  on  the  two  sides  of  the  boundary  would  be  discontinuous.)  The  only  other  quantity 
in  Eq.  23  that  is  directly  affected  by  the  boundary  is  y>[/Ax,  (/  -  l)Ay,  f\.  By  the  image 
method,  the  pressure  field  is  reflected  about  the  boundary,  so  that  y»[/Ax,  (/  -  l)Ay,  f\  = 
j3[/Ax,7Ay,  t\.  Therefore,  Eq.  23  becomes 


*The  term  “surface”  is  used  here,  although,  strictly  speaking,  the  rigid  condition  in  a  2-D  code  is 
implemented  along  a  line. 
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[(«  +  l/2)Ax:,7Aj,f]  +  V;,  [(f-l/2)Aj,7Aj,f]} 

X  {p[(/  +  \)^x,j^y,t]-p[{i-\)^x,j^y,t]] 

X  { [  /  Ax,  ( 7  +  1 )  Ay,  /]-/>[  /  Ax,  j Ay,  /  ] }  (60) 


-K  ( /Ax,  7  Ay,  / ) 


{ 


[ ( /  +  1  / 2 ) Ax, 7 Ay, /  ]-i4';f[(/-l/2) Ax, 7 Ay, / ] 


Ax 


[  /Ax,  ( 7  +  1  /  2 )  Ay,  /  ] 

A^^ 


} 


+/c(  /Ax,  7  Ay,  t)Q(  /  Ax,  7  Ay,  / ) . 


Similarly  Eq.  27  has  four  terms  involving  y  components  of  the  acoustic  particle 
velocity  and  wind  aty  =  (/  -  l/2)Ay,  all  of  which  are  zero.  Also  needed  are  the  x  compo¬ 
nents  of  the  acoustic  particle  velocity  and  wind  aty  =  (/  -  l)Ay.  But  because  the  x  coordi¬ 
nate  runs  parallel  to  the  boundary,  the  fields  for  these  quantities  are  also  reflected,  so 
Wx[(/  +  1/2) Ax,  (/  -  l)Ay,  /]  =  Wx[(/  +  1/2) Ax,  7 Ay,  /]  and  Vx[{i  +  1/2) Ax,  (/  -  l)Ay,  /]  =  Vx[{i 
+  l/2)Ax,7Ay,  /].  Finally,  in  Eq.  27  there  are  two  terms  involving  y  components  of  the 
velocities  aty  =  {j  -  l/2)Ay  that  must  be  set  to  zero. 


Sound  Propagation  in  a  Porous  Medium  and  Implementation 

Most  common  outdoor  ground  surfaces  cannot  be  modeled  satisfactorily  as  ideal, 
rigid  surfaces.  This  is  because  sound  energy  propagates  into  the  pores  of  the  ground, 
where  it  is  dissipated  by  viscosity  and  thermal  conduction.  Ground  surfaces  with  rela¬ 
tively  large,  open  pores,  such  as  snow,  absorb  much  of  the  sound  energy  incident  upon 
them.  Surfaces  with  small  pores,  such  as  cement  and  asphalt,  reflect  most  of  the  energy. 
The  acoustical  behavior  of  soils  is  intermediate  between  these  extremes.  For  complete¬ 
ness,  we  therefore  consider  in  this  report  FDTD  implementation  of  a  porous  ground  layer. 
We  show  that  the  ground  can  be  modeled  with  equations  very  similar  to  those  of  the  air. 

Morse  and  Ingard  (1968,  Eq.  6.2.22  and  6.2.23)  propose  the  following  phenomenol¬ 
ogical  equations  for  sound  propagation  in  a  porous  medium: 
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^  =  -(/c-^/Q)V-w,  (61) 


—  +  C7W  =  -Vp,  (62) 

where  a  is  the  static  flow  resistivity,  Kp  is  the  bulk  modulus  of  the  air  in  the  pores,  Pp  is 
the  density  of  air  in  the  pores,  and  O  is  the  porosity  (void  fraction)  of  the  medium.  The 
bulk  modulus  is  bounded  by  its  isothermal  value  pc^ly  and  its  adiabatic  value  p(?  (where  c 
indicates  the  adiabatic  sound  speed),  with  isothermal  conditions  prevailing  for  small 
pores  and/or  low  frequencies.  Subsequent  authors,  such  as  Attenborough  (1983),  have 
shown  that  when  the  properties  of  the  material  are  considered  in  bulk,  pp  has  a  compli¬ 
cated  frequency  dependence  and  in  fact  becomes  a  complex  quantity  in  the  frequency 
domain.  At  low  frequencies,  pp  may  be  replaced  by  where  is  called  the  structure 
constant  or  effective  density  factor.  Furthermore,  can  be  related  to  the  tortuosity  factor 
q,  which  describes  the  slanting  of  the  pores.  (For  pores  at  a  fixed  slant  angle  d,  q  is 
simply  1/cos  6.)  The  relationship  between  Sc  and  q  can  be  determined  for  some  idealized 
cases,  such  as  circular  cylindrical  pores,  for  which  Sc  =  {4l3)q^lD.  (Attenborough  1983). 
However,  the  relationship  between  Sc  and  q  depends,  in  general,  on  the  geometry  of  the 
pores. 

Restricting  our  attention  to  frequencies  low  enough  that  pp  can  be  considered  a  real 
constant  and  the  bulk  modulus  is  isothermal,  we  may  rewrite  Eq.  6 1  and  62  slightly  to 
produce  the  following  model  equations  for  the  porous  medium: 

^  =  -KpV  ■  w,  (63) 


5w 

A  -^  +  o-w  =  -Vp, 


(64) 


where 


Kp  ^  pc^  l]a 


(65) 


and 


Pe  ^Spp/D. 


(66) 


FDTD  Simulation  of  Sound  Propagation 


23 


are  the  effeetive  bulk  modulus  and  density  of  the  fluid  saturating  the  porous  material, 
respeetively.  This  equation  set  is  nearly  the  same  as  the  one  used  by  Salomons  et  al. 
(2002),  exeept  that  they  use  Ke=  pc^lQ.,  whieh  would  apply  to  adiabatic,  as  opposed  to 
isothermal,  conduction  in  the  pores. 

The  main  distinction  between  Eq.  63  and  64  and  those  for  the  nonmoving  fluid  is  the 
presence  of  the  term  involving  the  static  flow  resistivity,  a.  Physically  this  term  implies 
that  a  static  pressure  gradient  will  not  accelerate  the  fluid  to  an  infinite  speed;  rather,  the 
fluid  velocity  approaches  an  equilibrium  that  is  determined  by  the  viscous  drag  of  the 
porous  frame.  Since  O  is  always  less  than  one  and  5^  is  always  greater  than  1,  the  effec¬ 
tive  density  pe  is  always  greater  than  its  value  in  the  fluid.  The  phase  speed, 

=  clyf}^  ,  is  always  less  than  c,  implying  that  the  wavelength  is  shorter  in  the 
porous  ground  medium  than  in  the  air.  Table  1  shows  some  typical  parameter  values  for 
porous  media. 


Table  1.  Typical  values  of  the  static  flow  resistivity, 
porosity,  and  tortuosity  for  common  porous 
ground  surfaces. 


Material 

Flow  resistivity,  a 
(Pa-s  m  ^) 

Porosity,  Q 

Tortuosity,  q 

Asphalt 

3  X  10^ 

0.1 

3.2 

Grass 

2  X  10® 

0.5 

1.4 

Forest 

1  X  10® 

0.6 

1.3 

Sand 

5  X  lO'* 

0.35 

1.6 

Snow 

1  X  10® 

0.6 

1.7 

Finite-difference  implementation  of  Eq.  63  is  the  same  as  an  implementation  for  a 
nonmoving  medium.  Implementation  of  Eq.  64  is  slightly  different.  Eet  us  begin  by 
rewriting  it  in  component  form  as 


8t 


-  -s^w^ 


(67) 


dWy 

dt 


=  -S^Wy 


(68) 


where  Sg  =  aipe  and  bg  =  Mpg.  On  the  spatially  staggered  grid,  we  would  approximate  Eq. 
67  with  a  finite  difference  centered  atx  =  {i  +  l/2)Aj  andy  =j^y,  as  discussed  in  Section 
3.  Similarly  the  approximation  for  Eq.  68  would  be  centered  atx  =  /Ax  andy  =  (/  +  1/2) Ay. 
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Assuming  be  and  Se  are  defined  at  the  pressure  grid  nodes,  we  can  determine  these 
quantities  at  the  nodes  where  they  are  needed  using  Eq.  6  and  7.  The  pressure  gradients 
follow  from  Eq.  1 0  and  1 1  as  before.  The  velocities  in  the  terms  involving  Sg  are  known 
directly  at  the  approximation  centers.  We  see  that  FDTD  implementation  of  the  porous 
layer  is  trivial  once  the  basic  procedure  for  the  moving  medium  has  been  worked  out. 

Absorbing  Boundary  Condition 

Most  FDTD  codes  use  artificial  absorbing  layers  (called  absorbing  boundary  condi¬ 
tions,  or  ABCs)  around  the  edge  of  the  domain  to  avoid  spurious  numerical  reflections  at 
the  edges.  Popular  methods  include  the  Cerjan  attenuating  layer  (Cerjan  et  al.  1985)  and, 
more  recently,  the  perfectly  matched  layer  (Berenger  1994,  Yuan  et  al.  1997).  An  inter¬ 
esting  question  is  whether  a  porous  layer  of  the  type  discussed  in  the  previous  section  can 
be  used  as  an  ABC.  As  a  starting  point,  consider  the  reflection  coefficient  for  waves  in  a 
fluid  incident  upon  a  solid  material  (e.g.,  Kinsler  et  al.  1982): 


j^_^c  cos  9i  -  pc  cos  0^ 
Zg  cos  0i  +  pc  cos  9,  ’ 


(69) 


where  is  the  characteristic  impedance  for  waves  propagating  in  the  solid,  9i  is  the  angle 
of  incidence,  and  9i  is  the  angle  of  transmission.  The  angles  9i  and  9t  are  measured  rela¬ 
tive  to  the  surface  normal.  Note  that  9t  is,  in  general,  complex,  with  the  imaginary  part 
being  associated  with  evanescent  interface  waves.  If  the  reflection  coefficient  is  to  be 
zero,  we  must  therefore  have 


Zg  _  cos  9t 
pc  cos  9i 


(70) 


In  addition,  Snell’s  law  implies  that 


k  sin  9i  =  kg  sin  9^ , 


(71) 


where  k  =  cole  is  the  wavenumber  in  the  fluid  and  kg  is  the  complex  wavenumber  in  the 
solid.  Using  Snell’s  law  to  eliminate  cos0r  from  Eq.  70  results  in 


k 


cos2  0,+  ^  sin^^,-  =1. 

pej  Uc 


(72) 
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Expressions  for  the  characteristie  impedance  and  complex  wavenumber  can  be 
derived  directly  from  Eq.  63  and  64.  We  look  for  plane-wave  solutions  to  these  equations 
and  therefore  set 


p  =  ^  ^  ^  0.  (73) 

Such  a  solution  represents  waves  traveling  in  the  +x  direction.  By  substitution  into 
Eq.  63  and  64,  we  find 

^  _  ^Pc  (74) 

B  CD  kc  ' 

where  Pc  =  Pe~^  is  called  the  complex  density.  Solving  the  preceding  equations  for 
the  complex  wavenumber,  we  have 


By  definition,  =  AIB,  which  is  now  determined  as 


=  ^Jl^ePc-  (76) 

Substitution  into  Eq.  72  now  yields 


^^cos^  Bi  +  -^sin2  0,.  =  1. 

P  PcC 


(77) 


If  we  were  to  set  and  pc=  p  vcc  the  porous  layer,  then  Eq.  77  would  be  an 

identity  for  all  angles  of  incidence.  This  outcome  is  to  be  expected,  since,  if  the  air  layer 
and  porous  medium  have  identical  properties,  there  is  in  effect  no  interface.  More  inter¬ 
esting  is  what  happens  when  there  is  attenuation  in  the  porous  medium  and  pc  is  complex. 
Writing  out  the  real  and  imaginary  components  of  Eq.  77  explicitly,  we  have 


{p,+ialQ})K, 

2  2 
P 


2  /I  '^eiPe  -  .  2  ^ 

cos^  0:  +  ,  \  , - W^sm^  6:  =1. 

(pl+a^lco^y 


(78) 
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If  a/io  is  very  small  in  comparison  to  p^,  this  relationship  can  be  satisfied  approxi¬ 
mately  for  all  9i  by  setting  Ke  =  pc'  and pe  =  p,  or  equivalently  5c  =  O  =  y  =  1 .  Thus,  we 
may  use  a  porous  layer  as  an  ABC  so  long  as  the  artificial  static  flow  resistivity  is  small 
enough  to  satisfy  the  relationship  a/wp  «  1  throughout  the  frequency  range  of  interest. 
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6  IMPLEMENTATION  AND  EXAMPLE  CALCULATIONS 


Numerical  Issues 

Dependence  of  time  step  on  Mach  number 

For  numerical  stability  the  time  step  must  satisfy 

Ar 

At<— ,  (79) 

u 

where  u  is  the  speed  at  which  the  energy  propagates  and  Ar  the  grid  spacing.  On  a  two- 
dimensional  grid,  Ar  =  \/ .  In  the  FDTD  literature,  it  is  common  to 
define  the  Courant  number  as 


(80) 


Therefore,  to  satisfy  Eq.  79,  we  must  satisfy  the  so-called  Courant  condition,  C  <  1. 
The  time  step  and  grid  spacing  must  be  chosen  with  this  constraint  in  mind.  Since  the 
grid  spacing  should  generally  be  a  small  fraction  of  a  wavelength  for  good  numerical 
accuracy,  the  Courant  condition  in  practice  imposes  a  limitation  on  the  maximum  time 
step  possible  for  stable  calculations.  It  should  also  be  kept  in  mind  that  satisfaction  of  the 
Courant  condition  does  not  ensure  accurate  calculations. 

Let  us  consider  the  implications  of  the  Courant  condition  for  propagation  in  a 
uniform  flow.  In  this  case,  u  is  determined  by  a  combination  of  the  sound  speed  and  wind 
velocity.  In  the  downwind  direction,  we  have  m  =  m+  =  c  -l-  v.  In  the  upwind  direction,  u  = 
U-  =  c-v.  The  wavelengths  in  these  two  directions  are  l+  =  {c  +  v)//'and  I-  =  (c  -  v)lf, 
respectively,  where /is  the  frequency.  Since  the  wavelength  is  shortest  in  the  upwind 
direction,  it  is  that  dictates  the  grid  spacing.  We  set 


Ar  =  — =  A(i_m), 

N  ’ 


(81) 


where  N  is  the  number  of  grid  points  per  wavelength,  M  =  vie  is  the  Mach  number,  and  I 
=  Cz/is  the  wavelength  for  the  medium  at  rest.  We  see  that  a  finer  grid  is  required  as  M 
increases.  Regarding  the  time  step,  the  Courant  condition  implies 
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A,  <4.  (82) 

This  condition  is  most  difficult  to  meet  when  u  is  largest,  which  is  the  case  in  the 
downwind  direction.  Therefore,  we  must  use  u+  in  the  preceding  inequality  if  we  are  to 
have  accurate  results  throughout  the  domain;  specifically,  we  must  set 


/l_  _  1  1-M 


(83) 


From  this  equation  it  follows  that  the  time  step  at  M  =  1/3  must  be  1/2  the  value  nec¬ 
essary  at  M  =  0.  At  M  =  2/3,  the  time  step  must  be  1/5  the  value  at  M  =  0.  And  at  M  =  1 ,  a 
vanishing  time  step  is  required.  The  shortening  of  the  required  time  step  and  grid  spacing 
combine  to  make  calculations  at  large  Mach  numbers  more  computationally  intensive. 


Algorithm  discussion 

The  various  algorithms  in  Section  4  differ  in  their  memory  requirements  and  number 
of  calculations  required  per  time  step.  In  the  following  we  consider  the  computational 
aspects  of  each  algorithm  individually. 

Euler  method,  Eq.  29-31.  This  method  is  first  order  in  At  and  involves  one  calcula¬ 
tion  of  the  time-derivative  functions  fp,fx,  and^,  per  time  step.  Regarding  memory  usage, 
the  fields  themselves  must  be  stored,  and  also  the  calculations  of  the  time-derivative 
functions  at  each  grid  point.  Although  one  might  initially  think  that  the  fields  can  be 
updated  in  situ  while  calculating  the  time-derivative  functions  on  a  grid-point-by-grid- 
point  basis  (thereby  not  requiring  storage  of  at  every  grid  point),  it  must  be  remembered 
that  the  derivatives  depend  on  neighboring  field  points  as  well  as  the  FD  center  of  the 
approximation.  Therefore,  if  the  fields  at  nearby  grid  points  have  been  overwritten,  the 
calculations  of  fp,fx,  and  fy  will  be  incorrect.  We  characterize  the  calculation  burden  of  the 
Euler  method  as  one  “unif  ’  per  time  step,  where  one  calculation  unit  is  defined  as  the 
time  necessary  for  a  single  evaluation  of  the  time-derivative  functions.  The  memory 
usage  is  two  “units”  per  time  step,  where  a  single  memory  unit  is  defined  as  the  memory 
required  to  store  all  of  the  fields  at  each  grid  point. 

Unstaggered  leap-frog  scheme,  Eq.  33.  Since  it  is  based  on  centered  approxima¬ 
tions,  this  method  is  second  order  in  At.  It  involves  one  calculation  of  the  time-derivative 
functions  fp,fx,  and  fy  per  time  step.  This  method  is  readily  implemented  by  storing  the 
fields  over  two  time  steps  (/  and 7  -1)  and  the  time  derivative  functions  at  one  time  step 
(/).  Actually,  a  clever  algorithm  could  avoid  the  need  for  calculating  and  storing^,^,  and 
fy  at  each  grid  point,  since  the  fields  at  time  step  j  -  \  can  be  overwritten  without  affecting 
the  calculation  of fp,fx,  and  f,  at  time  step  j.  But  such  a  procedure  would  introduce,  in 
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practice,  more  complicated  memory  operations.  So  we  characterize  the  unstaggered  leap¬ 
frog  scheme  as  having  a  calculation  burden  of  one  unit  per  time  step  and  memory  usage 
of  two  to  three  units  per  time  step,  depending  on  the  implementation. 

Runge-Kutta  2nd  order,  Eq.  34-38.  As  indicated  by  the  name,  this  method  is  sec¬ 
ond  order  in  At.  Two  calculations  of fp,fx,  and^.  per  time  step  are  needed  per  time  step.  In 
addition  to  storing  the  fields  (which  can  be  updated  in  situ)  and  time-derivative  functions 
at  each  time  step,  one  must  store  the  time-derivative  functions  estimated  for  the  end  of 
the  time  step.  Therefore,  this  method  has  a  computational  burden  of  two  units  and  mem¬ 
ory  usage  of  three  units  per  time  step. 

Runge-Kutta  4th  order,  Eq.  39-43.  This  method  is  fourth  order  in  At.  Four  calcula¬ 
tions  of fp,fx,  and  fy  are  needed  per  time  step.  In  addition  to  storing  the  fields  (which  can 
be  updated  in  situ),  four  sets  of  time-derivative  functions  are  needed  per  time  step. 
Actually,  if  one  calculates  the  four  time-derivative  functions  sequentially  and  updates  the 
fields  after  each  update,  only  three  sets  of  arrays  encompassing  the  computational  grid 
need  to  be  stored.  Therefore,  this  method  has  a  computational  burden  of  four  units  and 
memory  usage  of  three  to  five  units  per  time  step,  depending  on  the  implementation. 

Staggered  leap-frog  scheme,  Eq.  51-53  (and  discussion  in  the  paragraph  after 
these  equations).  This  method  is  second  order  in  At  for  the  terms  not  involving  motion 
in  the  propagation  medium  and  first  order  in  At  for  the  terms  involving  motion.  It  is  the 
same  as  the  Euler  scheme  in  its  computational  and  memory  requirements. 

Aldridge  staggered  scheme,  Eq.  54-59.  This  method  uses  centered  differences  for 
all  terms  to  achieve  second-order  accuracy  in  At.  One  evaluation  of  the  time-derivative 
functions  is  required  for  each  whole  time  step.  The  method  is  most  readily  implemented 
by  storing  the  complete  fields  over  two  time  steps,  with  additional  memory  needed  to 
store  the  complete  time-derivative  functions  at  one  time  step.  A  reduction  in  memory 
usage  is  possible  by  observing  that  when  the  velocities  are  updated,  the  velocities  are 
only  needed  at  one  previous  time  step;  similarly,  when  the  pressure  is  updated,  the  pres¬ 
sure  is  needed  only  at  one  previous  time  step.  Therefore,  some  overwriting  is  possible. 
This  method  has  a  calculation  burden  of  one  unit  per  time  step  and  memory  usage  of  two 
to  three  units  per  time  step,  depending  on  the  implementation. 

The  results  of  this  discussion  are  summarized  in  Table  2.  In  comparing  the  methods, 
it  should  be  kept  in  mind  that  the  Euler  method  is  in  actuality  unstable  and  therefore 
cannot  really  be  considered.  The  staggered  leap-frog  scheme  is  also  of  questionable 
validity.  One  might  be  tempted  to  rule  out  the  Runge-Kutta  second-order  method,  since  it 
requires  twice  as  many  calculations  per  time  step  as  the  unstaggered  leap-frog  and  the 
Aldridge  staggered  methods.  However,  methods  of  the  same  temporal  order  do  not  neces¬ 
sarily  possess  the  same  accuracy.  If  one  can  take  time  steps  with  the  second-order  Runge- 
Kutta  that  are  twice  as  long  while  maintaining  the  same  accuracy  as  the  unstaggered  leap¬ 
frog  and  the  Aldridge  staggered  methods,  then  the  additional  computational  burden  per 
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time  step  would  be  worthwhile.  Therefore,  a  more  careful  study  of  the  properties  of  the 
methods  is  required.  We  will  consider  some  empirical  comparisons  between  the  various 
methods  later  in  this  section.  A  more  rigorous  theoretical  analysis  is  not  considered  in 
this  report,  however. 


Table  2.  Comparison  of  some  methods  for  advancing  in  time  the 
equations  of  sound  propagation  in  a  moving  medium.  One  unit  of 
memory  usage  is  defined  as  the  storage  required  for  arrays  of  all  the  field 
components  (acoustic  pressure  and  particle  velocities)  at  all  grid  nodes. 
One  unit  of  computational  burden  is  defined  as  the  number  of  operations 
required  to  evaluate  the  time-derivative  functions  once  for  all  of  the 
fields. 


Method 

Order  in  At 

Memory  usage 

Computational  burden 

Euler 

1 

2 

1 

Unstaggered  leap-frog 

2 

2-3 

1 

Runge-Kutta,  second-order 

2 

3 

2 

Runge-Kutta,  fourth-order 

4 

3-5 

4 

Staggered  leap-frog 

1-2 

2 

1 

Aldridge  staggered 

2 

2-3 

1 

Implementation  as  a  Matlab  code 

We  have  implemented  the  two-dimensional  propagation  equations  and  numerical 
methods  described  in  Sections  3-5  in  a  serial-processing  Matlab  code.  The  code  provides 
a  convenient  environment  for  testing  numerical  methods  for  FDTD  calculations  in  a 
moving  medium.*  The  time  derivatives  of  the  fields  (as  given  by  Eq.  23,  26,  and  27)  are 
implemented  explicitly  as  separate  functions,  allowing  all  of  the  temporal  finite- 
difference  schemes  discussed  in  Section  4  to  be  implemented  with  ease. 

Spatial  finite  differencing  is  accomplished  in  the  code  using  the  Matlab  circshift 
function,  which  circularly  shifts  numerical  arrays  upward,  downward,  rightward,  or  left¬ 
ward.  For  example,  a  circular  shift  downward  would  move  all  rows  down  one  position  in 
the  array,  except  for  the  last  row,  which  is  moved  to  the  top  of  the  array.  In  the  finite  dif¬ 
ference  code,  this  means  that  a  wave  propagating  downward  through  the  lower  domain 
boundary  reappears  at  the  upper  boundary.  This  non-physical  behavior  must  be  kept  in 
mind  when  determining  the  total  time  of  a  simulation,  so  that  spurious  “wrap-around” 


*A  full  3-D,  Fortran  code  that  implements  all  of  the  important  features  in  the  2-D  Matlab  code  is 
now  under  development.  The  3-D  code  performs  calculations  in  parallel  by  segmenting  the 
domain  into  rectangular-shaped  boxes  and  using  MPI  (Message  Passing  Interface)  instructions  to 
exchange  fields  at  common  grid  boundaries. 
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solutions  are  not  present.  Alternatively  an  absorbing  porous  layer  may  be  inserted  in  the 
vicinity  of  a  boundary  to  prevent  the  wrap-around. 

The  2-D  Matlab  code  has  also  been  written  to  accept  user-specified  perfectly  rigid 
and  acoustically  porous  objects  in  arbitrary  number  and  positions.  Technically  the  model 
for  rigid  surfaces  in  Section  5  assumed  a  surface  of  infinite  extent.  However,  the  code 
will  allow  users  to  specify  finite-sized  rigid  objects.  Corners  are  handled  by  applying  the 
infinite-surface  model  right  up  to  the  comer. 

No  Flow  with  Rigid  Boundary 

To  test  the  FDTD  code  and  rigid-boundary  implementation,  we  consider  in  this  sec¬ 
tion  a  set  of  calculations  for  a  harmonic  source  above  a  rigid  boundary  in  a  stationary 
medium.  Analytical  solutions  to  this  problem  are  known.  Specifically,  for  a  point  source 
radiating  into  an  infinite  2-D  medium  (or  equivalently  a  cylinder  radiating  into  a  3-D 
medium),  Morse  and  Ingard  (1968,  Sec.  7.3)  show 

f  +{y-ys  f  ),  (84) 

where  i  -  Vm  ,  (xs,ys)  is  the  source  position,  k=  Injlc  is  the  wavenumber, /is  the  source 
frequency,  and  is  the  zeroth-order  Hankel  function.  The  bolding  of  p  indicates  the 
complex  phasor  of  the  pressure  field,  as  opposed  to  its  real  value.  The  rigid  surface  can 
be  handled  by  the  method  of  images  (Morse  and  Ingard  1968,  Sec.  7.4),  yielding 


Pix,y)  =  i^k^]{x-x,  f  +(t-t,  )^  ) 

f  +{y  +  ys  f  )• 


(85) 


In  the  simulations,  we  implement  the  harmonic  source  with  a  finite  duration  signal 
that  is  tapered  at  the  initiation  and  termination.  The  tapering  alleviates  numerical  inaccu¬ 
racies  that  are  present  when  there  is  an  abmpt  change  in  the  source  emission.  A  cosine 
window  taper  is  used  here,  as  suggested  to  the  authors  by  N.  Symons.  The  tapered  source 
equation  is 


6(0- 


(^/2)[l  -  cos(;rf/ri  )]cos(2;r/  +  ^), 
Acos{27Tf  +  ^), 

(^/2)[1  -t  cos(;r(f  -  T)/T2  )]cos(2;r/  +  (/>), 


0, 


0  <  f  <  /, 
Ti<t<T-T2, 
T-T2<t<T, 
otherwise. 


(86) 
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In  the  preceding  equation,  A  is  the  source  amplitude,  ^^is  the  source  phase,  Ti  is  the 
duration  of  the  initiation  taper,  and  T2  is  the  duration  of  the  termination  taper.  All  calcu¬ 
lations  in  this  report  use  tapering  over  an  interval  of  three  periods  in  the  harmonic  wave 

{Ty  =  T2  =  2>IA 

The  domain  configuration  for  the  no-flow  calculations  is  shown  in  Fig.  2.  The  overall 
dimensions  of  the  domain  are  200  m  in  the  horizontal  and  1 00  m  in  the  vertical.  The 
origin  is  taken  to  be  the  middle  of  this  domain.  A  100-Hz  source  is  placed  at  (-60  m,  -30 
m).  A  rigid  surface  is  present  along  the  entire  lower  grid  boundary.  Absorbing  layers,  20 
m  thick,  are  present  at  the  top  and  left  sides  of  the  computational  grid.  These  absorbing 
layers  have  5^  =  O  =  y  =  1  as  suggested  earlier.  The  static  flow  resistivity  is  tapered  from 
100  Pa-s  m“^  at  the  side  of  the  layer  closest  to  the  source  to  1000  Pa-s  m“^  at  the  other 
side. 


Time  =  0.000  s 


Figure  2.  Configuration  of  the  computationai  domain  for  caicuiations  of  the  fieid 
generated  by  a  harmonic  source  above  a  rigid  boundary. 

The  calculations  in  this  section  use  the  second-order  spatial  differencing  scheme 
discussed  earlier  in  this  report,  together  with  a  fourth-order  Runge-Kutta  scheme  for 
marching  the  solution  forward  in  time.  Two  numerical  grid  resolutions  are  considered:  a 
low-resolution  grid  with  600  x  300  points  and  a  high-resolution  grid  with  1200  x  600 
points.  Since  the  wavelength  of  sound  at  100  Hz  is  3.4  m  in  air,  about  10  grid  points  per 
wavelength  have  been  used  in  the  low-resolution  grid  and  20  grid  points  per  wavelength 
in  the  high-resolution  grid.  The  Courant  number  (Eq.  80)  was  set  to  0.8,  leading  to  a  time 
step  of  0.58  ms  for  the  low-resolution  grid  and  0.29  ms  for  the  high-resolution  grid. 
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Figure  3  shows  the  pressure  field  0.450  s  after  initiation  of  the  100-Flz  source.  At  this 
time,  the  waves  have  propagated  approximately  150  m  from  the  source  in  the  +x  direc¬ 
tion.  Rapid  attenuation  of  the  sound  energy  is  seen  in  the  other  directions,  where  absorb¬ 
ing  layers  are  present.  A  series  of  spoke-like  pressure  minima  radiating  outward  from  the 
source  are  clearly  evident.  These  are  caused  by  destructive  interference  between  the 
direct  waves  and  the  waves  reflected  from  the  rigid  lower  boundary. 
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Figure  3.  Sound  pressure  field  0.450  s  after  initiation  of  the  100-Hz  sound  source. 

Comparisons  of  calculated  sound  transmission  loss  (TL)*  to  the  theoretical  result,  Eq. 
85,  are  presented  in  Fig.  4  and  5.  These  figures  show  results  for  the  low-  and  high- 
resolution  grids,  respectively.  The  receiver  in  this  case  is  at  the  same  height  as  the  source, 
and  the  TL  is  plotted  for  locations  increasing  in  horizontal  distance  from  the  source.  The 
“dips”  in  the  TL  result  from  the  interference  minima  mentioned  in  the  preceding  para¬ 
graph.  The  sound-pressure  amplitude  needed  for  the  TL  calculations  was  determined  by 
locating  the  maximum  absolute  value  of  the  pressure  within  the  final  30  time  steps  (one 
to  two  wave  periods)  of  the  received  signal.  Very  close  agreement  is  obtained  between 
the  theory  and  the  high-resolution  FDTD  calculations.  The  interference  minima  in  the 
low-resolution  calculations  are  shifted  to  noticeably  smaller  ranges  than  predicted  by  the 
theory. 

*Sound  transmission  loss  is  the  ratio  of  the  pressure  amplitude  recorded  at  the  receiver  to  the 
hypothetical  value  that  would  be  observed  in  free  space  at  a  distance  of  1  m  from  the  source.  The 
TL  is  expressed  in  decibels  by  taking  the  base-ten  logarithm  of  the  ratio  and  multiplying  by  20. 
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1.  Comparison  of  the  numerically  calculated  sound  transmission  loss  with 
tretical  result  for  the  low-resolution  numerical  grid. 
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Range  (m) 

Figure  5.  Comparison  of  the  numerically  calculated  sound  transmission  loss  with 
the  theoretical  result  for  the  high-resolution  numerical  grid. 


Propagation  over  Porous  Ground 

We  next  eonsider  FDTD  ealeulations  with  a  porous  ground  model  given  by  Eq.  63 
and  64.  The  domain  eonfiguration  (Fig.  6)  is  similar  to  the  one  used  for  the  no-flow 
ealeulations  (Fig.  2),  exeept  that  the  souree  has  been  moved  to  the  eenter  of  the  domain,  a 
rigid  barrier  has  been  plaeed  20  m  to  the  right  of  the  souree,  and  a  20-m-thiek  ground 
layer  is  now  present  at  the  bottom  of  the  domain.  The  absorbing  layer  to  the  left  of  the 
souree  has  been  removed.  The  souree  in  the  simulation  emits  1 0  eyeles  of  a  1 00-FIz  sine 
wave,  with  tapering  applied  during  the  three  eyeles  after  initiation  and  preeeding  termi¬ 
nation.  The  rigid  barrier  is  implemented  using  the  surfaee  boundary  eondition  deseribed 
for  rigid  surfaees  in  Seetion  5. 
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Time  =  0.000  s 


x(m) 

Figure  6.  Configuration  of  the  computationai  domain  for  caicuiations  invoiving 
porous  ground  and  a  rigid  barrier. 

Figures  7-9  show  calculations  for  an  acoustically  very  soft  surface.  The  porous  mate¬ 
rial  parameters  are  c  =  100  Pa-s  m“^,  q  =  1.8,  and  O  =  0.5.  In  the  first  figure  of  the 
sequence,  taken  0.057  s  after  source  initiation,  the  sound  waves  are  just  beginning  to 
impinge  on  the  barrier.  By  the  second  figure  (Fig.  8),  taken  at  0. 154  s,  a  wavetrain  is  seen 
reflecting  strongly  off  the  barrier  while  a  second,  weaker  wavetrain  reflects  from  the 
ground.  A  diffracted  wavetrain  curls  around  the  barrier.  The  shorter  wavelength  of  the 
sound  and  strong  attenuation  within  the  ground  are  readily  evident.  In  Figure  9,  taken  at 
0.250  s,  the  initial  wave  and  barrier  reflection  are  clearly  seen  propagating  leftward.  The 
wavetrain  that  propagated  over  the  top  of  the  barrier  (including  the  diffracted  wave) 
propagates  to  the  right.  Weak  ground  reflections  associated  with  each  of  these  three 
wavetrains  are  also  visible.  Waves  in  the  ground  refract  strongly  toward  the  surface 
normal,  as  is  expected  from  Snell’s  law. 
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Figure  7.  Propagation  above  soft,  porous  ground  in  the  presence  of  a  rigid  barrier, 
0.057  s  after  source  initiation. 
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Figure  8.  Propagation  above  soft,  porous  ground  in  the  presence  of  a  rigid  barrier, 
0.154  s  after  source  initiation. 
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Time  =  0.249  s 


x(m) 

Figure  9.  Propagation  above  soft,  porous  ground  in  the  presence  of  a  rigid  barrier, 
0.249  s  after  source  initiation. 

By  comparison,  a  calculation  for  a  moderately  hard  ground  surface  is  shown  in 
Figure  10.  This  calculation  is  identical  to  the  preceding  one,  except  that  a  has  been 
increased  to  1 0“"^  Pa-s  m"^.  This  value  is  characteristic  of  very  soft  soil.  The  pressure  field 
0.250  s  after  source  initiation  is  shown.  Notice  that  the  higher  value  for  the  static  flow 
resistivity  causes  the  sound  waves  to  reflect  nearly  perfectly  from  the  ground  in  this  case. 
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Time  =  0.250  s 


x(m) 

Figure  10.  Propagation  above  moderateiy  hard,  porous  ground  in  the  presence  of  a 
rigid  barrier,  0.250  s  after  source  initiation. 

Calculations  involving  porous  media  with  flow  resistivities  exceeding  1 0“"^  were  also 
attempted  but  became  unstable  when  the  sound  wave  impinged  on  the  porous  layer.  Since 
the  porous  layer  is  highly  attenuative  for  large  flow  resistivities,  the  instability  could  have 
been  caused  by  poor  resolution  of  the  spatial  scales  over  which  the  wave  attenuates. 
Further  research  will  be  necessary  to  better  understand  the  grid  spacings  and  time  steps 
necessary  to  achieve  stability  in  a  highly  dissipative  porous  medium. 

High  Mach  Number,  Uniform  Flow 

In  this  section,  calculations  are  considered  for  very-high-speed  uniform  flows. 
Although  the  flow  speeds  are  much  higher  than  would  normally  be  encountered  in  the 
atmosphere,  the  calculations  are  very  useful  for  testing  the  finite-difference  implementa¬ 
tion  of  the  terms  in  Eq.  1  and  2  that  are  particular  to  a  moving  medium.  The  source  is  the 
same  tapered,  1 0-cycle,  1 00-Hz  signal  considered  in  the  previous  section.  The  domain 
size  for  all  of  the  calculations  is  1 00  x  1 00  m.  Most  of  the  calculations  have  800  grid 
points  in  each  direction,  although  a  higher-resolution  calculation  with  1200  grid  points  is 
also  discussed.  Calculations  at  three  Mach  numbers  are  considered:  0,  0.3  (about  104 
m/s),  and  0.6  (about  207  m/s).  The  number  of  grid  points  per  wavelength,  in  the  upwind 
direction,  is  27.2,  19.0,  and  10.9,  respectively,  for  these  three  Mach  numbers.  The  fourth- 
order  Runge-Kutta  scheme  with  a  Courant  number  of  0.4  was  used  to  advance  the  wave- 
fields  in  time. 
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Figures  1 1-13  show  the  full,  2-D  pressure  fields  at  Maeh  0,  0.3,  and  0.6  after  0. 1 1  s 
(0.01  s  after  the  souree  has  terminated).  Elongation  of  the  waves  propagating  in  the 
direetion  of  the  flow,  and  eompression  of  the  waves  propagating  eounter  to  the  flow,  are 
elearly  observed.  (In  the  Maeh  0.6  ealculation,  wavefronts  appear  on  the  left  edge  of  the 
figure.  These  are  wrapped  around  from  the  right  edge  by  the  periodie  BC.) 
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Figure  11.  Sound  pressure  field  in  the  absence  of  a  back¬ 
ground  flow,  Mach  0  flow. 
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Mach  0.3,  Time  =  0.110  s 
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Figure  12.  Sound  pressure  field  for  propagation  in  a  uni¬ 
form,  Mach  0.3  flow.  The  flow  direction  is  left  to  right. 
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Figure  13.  Sound  pressure  field  for  propagation  in  a  uni¬ 
form,  Mach  0.6  flow.  The  flow  direction  is  left  to  right. 

Exact,  time-domain  solutions  for  2-D  wave  propagation  in  a  moving  medium  are  not 
generally  available.  Recently,  V.E.  Ostashev*  has  calculated  the  dependence  of  the  pres¬ 
sure  amplitude  on  azimuth  for  a  monofrequency  sound  source  radiating  in  two  dimen¬ 
sions.  The  result  is 


p{r,e,M)  sin^  ef  -  M  COS0 

p{r,0,M  =  O)  0)^'\\-M) 


(87) 


where  p{r,  0,  M)  is  the  pressure  amplitude  as  a  function  of  cylindrical  distance  r,  azi¬ 
muthal  angle  0,  and  Mach  number.  This  equation  is  intended  for  application  in  the  far- 
field,  that  is,  many  wavelengths  from  the  source.  Figure  14  compares  the  FDTD 
calculations  in  the  previous  three  figures  to  Ostashev’s  theory.  Very  close  agreement  is 


Personal  communication,  2003. 
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obtained  at  Mach  0.3.  At  Mach  0.6  the  amplitude  from  the  FDTD  calculations  is  some¬ 
what  lower  than  expected  in  the  upwind  direction.  A  higher-resolution  grid  noticeably 
improves  the  predictions.  Taken  together,  these  results  demonstrate  that  excellent  agree¬ 
ment  with  Ostashev’s  theory  can  be  obtained,  although  a  very-high-resolution  grid  is 
required  for  high  Mach  numbers. 


Azimuth,  a  (deg) 

Figure  14.  Comparison  of  the  angular  dependence  of  the  pressure  amplitude  from 
the  FDTD  calculations  of  Ostashev’s  theory.  The  upwind  direction  is  0°  and  the 
downwind  direction  is  180°. 

Figures  15  and  16  repeat  the  Mach  0.3  calculation  shown  in  Figure  12  except  that 
different  temporal  integration  methods  are  used.  Results  for  the  “staggered  leap-frog” 
method  (Fig.  15)  are  qualitatively  similar  to  the  fourth-order  Runge-Kutta  results.  The 
time  step  used  for  the  staggered  leap-frog  calculations  was  1/4  (Courant  number  C  =  0.1) 
that  used  for  the  fourth-order  Runge-Kutta  calculations,  approximately  equalizing  the 
amount  of  computational  effort  involved  (Table  2).  Close  examination  of  the  near-source 
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wavefronts  propagating  in  the  downwind  direetion  reveals  a  noisy  appearance  not  present 
in  Figure  12.  This  noise  is  actually  a  growing  numerical  instability  in  the  calculations. 
Such  numerical  instability  is  even  more  pronounced  in  the  second-order  Runge-Kutta 
calculations  (Fig.  16).  These  calculations  were  performed  with  C  =  0.2. 
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Figure  15.  Sound  pressure  field  for  propagation  in  a  uni¬ 
form,  Mach  0.3  flow.  The  staggered  leap-frog  method  was 
used  to  advance  the  solution  in  time. 
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Mach  0.3,  Time  =  0.110  s 


x(m) 


Figure  16.  Sound  pressure  field  for  propagation  in  a  uni¬ 
form,  Mach  0.3  flow.  The  second-order  Runge-Kutta  method 
was  used  to  advance  the  solution  in  time. 

Downwind  and  upwind  results  at  Mach  0.3  are  compared  for  several  temporal  inte¬ 
gration  methods  in  Figures  17  and  18.  Shown  are  fourth-order  Runge-Kutta  (C  =  0.4), 
unstaggered  leap-frog  (C  =  0.1),  and  staggered  leap-frog  (C  =  0.1).  The  first  two  of  these 
methods  yield  graphically  indistinguishable  results.  The  staggered  leap-frog  scheme, 
however,  yields  a  signal  with  a  noticeably  larger  amplitude.  The  numerical  instabilities 
are  also  clearly  evident  at  small  distances  from  the  source.  The  lack  of  agreement 
between  the  staggered  leap-frog  scheme  and  the  other  methods,  combined  with  its  non¬ 
rigorous  method  of  derivation,  strongly  suggests  that  it  is  inferior  for  this  problem. 
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Downwind 


Figure  17.  Sound  pressure  traces  for  downwind  propagation  at  Mach  0.3.  Shown 
are  resuits  for  the  fourth-order  Runge-Kutta,  staggered  ieap-frog,  and  unstaggered 
ieap-frog  methods. 
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Upwind 


Figure  18.  Sound  pressure  traces  for  upwind  propagation  at  Mach  0.3.  Shown  are 
results  for  the  fourth-order  Runge-Kutta,  staggered  leap-frog,  and  unstaggered 
leap-frog  methods. 

A  similar  comparison  of  downwind  and  upwind  results,  except  at  Mach  0.6,  is  shown 
in  Figures  19  and  20.  The  temporal  integration  methods  shown  in  these  figures  are  fourth- 
order  Runge-Kutta  (C  =  0.4),  unstaggered  leap-frog  (C  =  0.1),  and  Aldridge  (C  =  0.1). 

The  second-order  Runge-Kutta  and  staggered  leap-frog  methods  were  both  extremely 
unstable  at  this  Mach  number  and  therefore  are  not  shown.  The  three  methods  displayed 
in  the  figure,  however,  were  all  stable  and  provide  graphically  indistinguishable  results. 


Sound  pressure  (Pa) 
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Downwind 


Distance  from  source  (m) 

Figure  19.  Sound  pressure  traces  for  downwind  propagation  at  Mach  0.6.  Shown 
are  results  for  the  fourth-order  Runge-Kutta,  Aldridge,  and  unstaggered  leap-frog 
methods. 
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Upwind 


Figure  20.  Sound  pressure  traces  for  upwind  propagation  at  Mach  0.6.  Shown  are 
resuits  for  the  fourth-order  Runge-Kutta,  Aidridge,  and  unstaggered  ieap-frog 
methods. 

Based  on  these  eomparisons  of  the  temporal  integration  methods,  it  appears  that  the 
unstaggered  leap-frog,  Aldridge,  and  fourth-order  Runge-Kutta  methods  all  provide  good 
aeeuraey  and  stability  for  a  similar  level  of  eomputational  effort.  More  detailed  analysis 
would  be  neeessary  to  determine  whieh  of  these  methods  is  most  satisfaetory.  From  a 
praetieal  standpoint,  any  of  them  appears  to  be  a  reasonable  ehoiee. 

Nonuniform  (Turbuient)  Fiow 

Atmospherie  turbulenee  has  several  important  effeets  on  signals  reeeived  by  aeoustie 
sensors:  it  eauses  signal  levels  to  fluetuate  randomly,  it  distorts  wavefronts  (thereby 
ehanging  angles-of-arrival  and  hindering  souree  direetion  finding),  and  it  seatters  sound 
energy  over  and  around  objeets.  Therefore,  incorporation  of  turbulence  into  FDTD 
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calculations  is  of  much  practical  importance.  Two  general  approaches  to  generating  the 
necessary  turbulence  input  fields  are  possible:  kinematic  and  dynamic.  The  kinematic 
method  consists  of  synthesizing  random  fields  that  possess  certain  desired  statistical 
properties  characteristic  of  actual  turbulence.  The  dynamic  method  consists  of  simulating 
actual  turbulence  with  fluid  dynamical  equations. 

The  primary  advantages  of  the  kinematic  methods  are  low  computational  burden,  the 
ability  to  readily  generate  fields  at  a  high  spatial  resolution,  and  tight  control  over  statisti¬ 
cal  properties.  This  last  feature  is  useful,  for  example,  in  comparing  the  FDTD  calcula¬ 
tions  to  known  analytical  solutions.  The  primary  advantage  of  the  dynamic  methods  is 
that  the  temporal  properties  and  the  higher-order  statistics  of  the  fields  are  as  close  as 
possible  to  the  actual  atmosphere.  The  main  dynamic  method  of  interest  for  FDTD  cal¬ 
culations  at  Army  tactical  scales  is  large-eddy  simulation,  or  LES  (Moeng  1998).  At  pre¬ 
sent,  LES  is  capable  of  providing  inputs  to  an  FDTD  code  at  a  resolution  of  several 
meters  (e.g.,  Sullivan  et  al.  1996).  Conceivably  the  EES  could  be  directly  coupled  with 
the  acoustical  calculation  by  simultaneously  solving  a  single  set  of  fluid  equations.  FIow- 
ever,  such  a  procedure  can  lead  to  poor  results  because  of  the  discrepancy  between  the 
dominant  scales  of  the  atmospheric  turbulence  process  and  the  acoustic  wavelength,  as 
well  as  the  small  magnitude  of  the  acoustic  fluctuations  relative  to  the  atmospheric  ones. 

For  illustrative  purposes  we  consider  here  a  single  example  of  FDTD  calculation  of 
sound  propagation  through  turbulence.  The  example  uses  the  recently  developed  tech¬ 
nique  for  generating  kinematic  turbulence  called  quasi-wavelets  (QWs).  Eike  ordinary 
wavelets,  QWs  are  based  on  rescaling  and  translation  of  a  specified  parent  function. 
Unlike  ordinary  wavelets,  the  positions  and  orientations  of  the  QWs  are  random. 
Goedecke  et  al.  (2003)  derived  a  random  QW  distribution  and  parent  function  that 
exactly  reproduce  a  prescribed  turbulence  spectrum  such  as  von  Karman’s.  Figures  21 
and  22  show  snapshots  of  and  Vy,  respectively,  for  a  random  QW  field.  As  is  charac¬ 
teristic  of  atmospheric  turbulence,  there  are  strong  but  infrequent  patches  of  large-scale 
velocity  fluctuations,  interspersed  with  more  numerous  but  weaker  small-scale  patches. 
FDTD  calculations  of  the  propagation  of  a  six-cycle  burst  at  1 00  FIz  through  the  QW 
fields  is  shown  in  Figure  23.  Note  the  increasing  distortions  to  the  wavefronts  as  they 
propagate  from  the  source.  This  example  is  somewhat  exaggerated  in  that  the  turbulent 
velocity  fluctuations  have  a  standard  deviation  of  about  10  m  s“’,  which  is  an  order  of 
magnitude  higher  than  would  typically  be  encountered  in  practice.  Flowever,  the  propa¬ 
gation  distance  is  much  shorter  than  typical  in  Army  tactical  scenarios. 
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Figure  21.  Snapshot  of  a  random  turbulence  field  generated  by  the  quasi¬ 
wavelet  method.  The  x  component  of  the  wind  velocity  is  shown. 
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Figure  22.  Same  as  Figure  21,  except  that  the  y  component  of  the  wind 
veiocity  is  shown. 
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Figure  23.  Propagation  of  a  six-cycle  burst  at  100  Hz  through  the  turbu¬ 
lence  fields  shown  in  Figures  21  and  22. 


54 


ERDC/CRRELTR-04-12 


7  CONCLUSION 

This  report  has  addressed  many  aspects  of  finite-difference,  time-domain  (FDTD) 
simulation  of  sound  wave  propagation  in  a  moving  medium,  from  the  theoretical  equa¬ 
tions  upon  which  it  is  based  to  practical  issues  of  numerical  implementation.  One  of  the 
most  important  insights  is  that  motion  in  the  propagation  medium  makes  impossible 
direct  implementation  of  the  customary,  staggered  grid  schemes  used  for  other  types  of 
wave  propagation  problems.  This  is  because  of  the  presence  of  the  advective  terms  in  the 
differential  equations,  which  do  not  naturally  center  themselves  on  the  grid  points  where 
they  are  needed.  Therefore,  alternative  finite- difference  formulations  must  be  considered. 

The  numerical  approaches  considered  in  this  report  were  all  based  on  the  customary, 
staggered  spatial  grid,  in  which  the  pressure  and  particle  velocities  are  displaced  by  one- 
half  of  the  grid  interval.  Advancement  of  the  wavefield  variables  in  time  was  then  con¬ 
sidered  in  terms  of  unstaggered  and  staggered  temporal  grids.  Two  temporally  unstag¬ 
gered  schemes — ^the  fourth-order  Runge-Kutta  method  and  the  unstaggered  “leap-frog” 
method — were  found  to  be  quite  satisfactory.  A  temporally  staggered  scheme,  with  fields 
stored  from  two  previous  time  levels,  also  performed  well.  However,  a  non-rigorous, 
staggered,  leap-frog  scheme  analogous  to  the  one  currently  in  widespread  use  for  FDTD 
wave  propagation  calculations  in  a  non-moving  medium  was  found  to  be  inaccurate. 
Additional  comparisons  of  the  various  temporal  integration  schemes  would  be  very 
worthwhile. 

Rigid  boundary  conditions  and  porous  ground  layers  were  also  considered  in  this 
report.  The  rigid  condition  is  useful  for  modeling  buildings,  barriers,  tree  trunks,  and 
other  acoustically  hard  objects.  A  simple,  low-frequency  model  for  the  porous  ground 
can  be  easily  incorporated  into  an  acoustic  FDTD  code.  However,  the  utility  of  this 
model  is  limited  to  materials  with  low  flow  resistivities,  such  as  snow  or  coarse  sand. 
Therefore,  it  would  be  highly  desirable  to  develop  other  approaches  to  implementing  the 
ground  material,  such  as  a  ground-layer  model  with  applicability  to  a  broader  frequency 
range  or  an  impedance  boundary  condition  that  does  not  require  an  actual  ground  layer. 
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